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We examine collective modes, stability, and BCS pairing in a quasi-two-dimensional gas of dipolar 
fermions aligned by an external field. By using the (conserving) Hartree-Fock approximation, which 
treats direct and exchange interactions on an equal footing, we obtain the spectrum of single-particle 
excitations and long wavelength collective modes (zero sound) in the normal phase. It appears that 
exchange interactions result in strong damping of zero sound when the tilting angle between the 
dipoles and the normal to the plane of confinement is below some critical value. In particular, zero 
sound cannot propagate if the dipoles are perpendicular to the plane of confinement. At intermediate 
coupling we find unstable modes that can lead either to collapse of the system or the formation of 
a density wave. The BCS transition to a superfluid phase, on the other hand, occurs at arbitrarily 
weak strengths of the dipole-dipole interaction, provided the tilting angle exceeds a critical value. 
We determine the critical temperature of the transition taking into account many-body effects as 
well as virtual transitions to higher excited states in the confining potential, and discuss prospects 
of experimental observations. 
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I. INTRODUCTION 

In recent years, experimentalists have achieved ma- 
jor breakthroughs in preparing samples of diatomic 
molecules in the rovibrational ground state and cooling 
them towards quantum degeneracy [fj-y] . With heteronu- 
clear molecules, in particular, rotational degrees of free- 
dom can be excited in a controlled way by applying ex- 
ternal electric fields, and are associated with large elec- 
tric dipole moments 0-0] ■ This possibility of inducing 
strong and anisotropic dipole-dipole interactions between 
the molecules opens fascinating prospects for the obser- 
vation of various many-body effects and novel quantum 
phases 0,11^. 

The above-mentioned experimental studies of het- 
eronuclear polar molecules P, [|| suffer from losses due 
to chemical reactions such as KRb + KRb — » K2 + Rb2 
[T7L [l8l ] , which place severe limitations on the achievable 
densities in three-dimensional samples. These reactions 
are significantly suppressed if one confines the molecules 
to a quasi-two-dimensional (quasi-2D) geometry and ori- 
ents their dipole moments perpendicular to the plane of 
the 2D translational motion |19l - [2l| , as has been verified 
experimentally by de Miranda et al. [22]. Moreover, in 
some of the polar molecules that consist of alkali atoms, 
atom-exchanging reactions are endothermic and, there- 
fore, do not occur [23]. Thus, it seems most promising 
for future investigations of dipolar molecules to either fo- 
cus on species that do not undergo chemical reactions, or 
consider samples that are strongly confined to a quasi-2D 
regime. 

In this paper, we consider a quasi-2D gas of fermionic 
dipoles, aligned by an external field (see Fig. [T]). The 
simplest case corresponds to the dipoles being oriented 



perpendicular to the plane of confinement, so that the 
pairwise dipole-dipole interaction is isotropic and repul- 
sive. At non-zero values of the tilting angle Oo it becomes 
anisotropic and, for 6q > arcsin(l/v / 3), partially attrac- 
tive. A discussion of this physical setup within the frame- 
work of Fermi liquid theory was given in [24j |. where, e. g., 
single-particle properties such as the anisotropic self- 
energy and the resulting deformation of the Fermi sur- 
face from the spherical shape corresponding to the non- 
interacting case were calculated to first order in pertur- 
bation theory. While this perturbative approach provides 
reliable answers in the weak coupling regime, moderate 
interaction strengths require more sophisticated methods 
such as the Hartree-Fock approximation (HFA), which 
was used in Ref. (25| to obtain the spectrum of single- 
particle excitations at zero and finite temperature. At 
zero temperature, the results were found to agree very 
well with the outcome of a variational approach that was 
initially used to study Fermi surface deformations in the 
3D cas e 12611 and adapted to the 2D case by the authors 
of Ref. pjj? 

Collective modes in single-, bi- and multilayered struc- 
tures of dipolar Fermi gases were studied in (27] using 
the random phase approximation (RPA), which neglects 
exchange interactions. As a result, in particular for the 
single-layer setup, these studies predict the spectrum of 
long wavelength collective excitations to be sensitive to 
microscopic details of the two-body interaction potential 
in the form of a short-distance cut-off, which is needed 
to handle the singular behavior of a dipole-dipole inter- 
action. In a quasi-2D setting, the characteristic length 
of the harmonic confinement takes the role of the cut-off, 
resulting in a confinement-dependent value of the RPA 
speed of zero sound [28|, H^] . The authors of [3Q] , how- 
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ever, correctly remark that a cut-off dependent constant 
term in the momentum space representation of the inter- 
particle potential corresponds to a short-range contact 
interaction in real space and, therefore, must not have 
an effect in a single-component Fermi gas. Thus, the 
existence of the RPA zero sound mode is questionable. 
In this paper, we study zero sound collective modes on 
the basis of the so-called conserving HFA developed in 
Refs. [3ll, HI]. The advantage of this method is that it 
provides a way to fully include exchange contributions in 
a given order of perturbation theory, such that the results 
are consistent with conservation of particle number, en- 
ergy, and momentum, as well as with fermionic statistics 
of particles. We show that the existence of a zero sound 
collective mode in a quasi-2D dipolar Fermi gas for small 
values of the tilting angle 80 and the dependence of the 
sound velocity on a short-distance cut-off are artifacts 
of the RPA. In particular, we find that the propagation 
of zero sound is not possible if the dipoles are aligned 
perpendicular to the plane of confinement or if they are 
tilted only slightly - in consistence with the homogeneous 
3D setting [2J, |33[ i n which there is no propagating zero 
sound in the directions perpendicular (or close to per- 
pendicular) to the direction of dipole polarization. 

The issue of stability of the normal phase of the system 
against collapse was addressed by Chan et al. 24| follow- 
ing Pomeranchuk's approach [3J|: The normal phase is 
thermodynamically stable if an arbitrary distortion of the 
Fermi surface results in an increase of the ground state 
energy. In Ref. 24], however, distortions do not refer to 
the deformed Fermi surface, but rather to the circular one 
of the non-interacting system, i. e., the authors of this ref- 
erence are performing the stability analysis around a con- 
figuration that does not extremize the ground state en- 
ergy. Hence their expression for the change in the ground 
state energy contains a term that is linear in the distor- 
tion [Eq. (81) in [24j], which is absent if one takes the 
deformed Fermi surface as reference center [see Eq. (jit?)) 
below]. An alternative approach by Bruun and Taylor 
[I2} uses a variational ansatz for the shape of the Fermi 
surface, on the basis of which the compressibility is cal- 
culated. The collapse instability is then identified with 
a negative value of this quantity. The stability of the 
normal phase against density fluctuations with a finite 
momentum (density wave instability) was investigated in 
RPA in Refs. [3(| and - extending the discussion to finite 
temperatures and taking into account the deformation of 
the Fermi surface - [25] . It was found, that a density wave 
transition takes place in a broad region in the parameter 
space (the coupling strength and the tilting angle 8 ), 
where the system is stable against collapse. An improve- 
ment of the RPA treatment was achieved in Ref. [35[ 
by an approximate treatment of correlations beyond the 
RPA in the form of a local-field correction in the density- 
density response function (the same method with a dif- 
ferent form of the local-field correction was used in (3(| 
for a specific value of the tilting angle 9q = arccos l/v3). 
The result is that the density wave instability should be 



expected at higher values of the interaction strength than 
predicted by the RPA. This approach, however, provides 
results that are very sensitive to the form of the local-field 
correction, and the justification of a particular choice is 
not clear. We address the issue of stability on the basis 
of the conserving HFA that provides a consistent (and 
within this approximation scheme exact) treatment of 
exchange contributions. For the case 9q = 0, the same 
approach was used in Ref. [13] to study the density wave 
instability in a dipolar mono- and multilayer sytems. Our 
result (see Sec. IVIIB[) for the critical value of the cou- 
pling strength for Qq = agrees very well with that from 
Ref. [13. 

The critical temperature T c of the transition to the 
superfluid phase for the quasi-2D dipolar Fermi gas was 
obtained by Bruun and Taylor [12j in the BCS approach 
with the dipole-dipole interaction restricted to the dom- 
inant p-wave channel. We extend this work by taking 
into account the full angular dependence of the dipole- 
dipole interaction, as well as by calculating the preex- 
ponential factor in the expression for the critical tem- 
perature. The latter requires taking into account both 
the many-body contributions to the interparticle inter- 
action (the so-called Gor'kov-Melik-Barkhudarov (GM) 
corrections 38]) and virtual transitions to excited states 
in the trapping potential, which ultimately result in a 
non-trivial dependence of the critical temperature on the 
trapping frequency, the gas density, and the tilting angle. 

This paper is organized as follows: In Sec.[H]we present 
the microscopic model for the quasi-2D dipolar Fermi 
gas and identify the relevant parameters and parame- 
ter regimes. We review key quantities in the many-body 
problem and the equations that relate them and discuss 
our strategy to solve these equations in the intermediate 
coupling regime in Sees. IIIII and IIV1 An investigation 
of the single-particle excitation spectrum in the normal 
phase is given in Sec. El and is a prerequisite for the 
study of collective modes in Sec. I VII We address the is- 
sue of instability towards, respectively, collapse and the 
formation of a density wave in Sees . I VIL^l and IVII B I Re- 
sults for the critical temperature of the transition to the 
superfluid phase are given in Sec. IVIIII Finally, Sec. IIXI 
is devoted to a summary of our findings and a discus- 
sion of the prospects of observing the described phenom- 
ena in experiments. Details of our numerical methods 
and analytical expressions for the matrix elements of the 
dipole-dipole interaction are given in the appendices. 



II. SYSTEM 

We consider a gas of dipolar fermions of mass m with 
dipole moments d = dd, which are polarized along the di- 
rection d = (sin 0o,O, cos 9q), i.e., 9o is the angle between 
the orientation of dipoles and the z-axis (see Fig.[T]). The 
gas is strongly confined to the xy-plane by a harmonic 
trapping potential V(z) — mwgZ 2 /2. Here, strong con- 
finement means that the transverse extension of the gas 
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As a result, first order interaction effects (see Sees. fVl 



FIG. 1. (Color online) Fermionic dipoles confined to the xy- 
plane. The dipoles are aligned in the rrz-plane and form an 
angle 9q with the z-axis. 



cloud, which is on the order of Iq = y/K/moJo, is small as 
compared to the mean interparticle separation in the xy- 
plane. The latter quantity is proportional to the inverse 
Fermi momentum Pp , which, in turn, is determined by 

the area density rc.223, Pf = hs/^Kn^D. Thus, the small 
parameter characterizing the tight-confinement or quasi- 
2D limit is r\ = Pp Iq/K <C I. This condition implies that 
the Fermi energy Bp is much smaller than the trapping 
potential level spacing, e F , /hu>o <C 1. 

We also assume ultracold temperatures, T <C Ep, such 
that the average kinetic energy of particles is given by the 
Fermi energy £p . 

The Hamiltonian of this system reads 



fi 2 
2m 



A H — muiiz 2 — u! 



+ - / drdr' ft(r)ft(r')V d (r - r')$(r')$(r), (I) 



where ip( r ) is the fermionic field operator, fj! denotes the 
chemical potential, and Vd(r) = (d 2 /r 3 )[l — 3(f-d) 2 ] is the 
dipole-dipole interaction. Here we omit the contribution 
of the short-range part of the interparticle interaction: In 
the considered case of a single-component Fermi gas, it 
results only in p-wave scattering, which is small assuming 

that p^ro/h <C 1, where r$ is the radius of the short- 
range part of the interparticle interaction. 

A characteristic length of the dipole-dipole interaction 
is given by = md 2 /h 2 3> tq. This is the length scale 
below which the dipole-dipole interaction substantially 
influences the relative wave function of two particles. We 
assume rd ^ la, such that interparticle collisions are es- 
sentially three-dimensional. This gives us another small 
parameter g = PpTd/h <C 1. 

Under the above conditions, the motion of particles 
in the z-direction is limited to the ground state of the 
confining potential <j>o(z) - the lowest harmonic oscil- 
lator (HO) level, and a single-particle wave function is 
ip(r) — <p(p) cf>o(z), where <p(p) describes the in-plane 
motion [p — (x,y)\. 



IVI1 and IVIip can be described by an effective interaction 

V Q (p) = J dzdz'Mz) 2 Vd(p,z-z')Mz') 2 (2) 
with the Fourier transform [p = (p x ,p y )} [28] 

V (p) = V^C^vfe) Hp)sin 2 o -2P 2 (cos0 o )] , 

where u(p) = (p 2 — Py)/p 2 and w(x) = xe x2 erfc(x)0 
When p s» Pp, the argument of w (^fJ is as small as 
plo /H~ rj <C 1 and we have 



Vb(p) « [ U (p)sin 2 o - 2P 2 (cos0 o ) 



(4) 



We note that in this limit the effective interaction is in- 
dependent of the confinement length Iq, 

Processes of second order in the interaction (see 
Sec. IVIlij) . on the other hand, involve virtual transitions 
to excited states <f> n (z) with n > of the harmonic con- 
fining potential. The matrix elements of the interaction 
for these transitions, in momentum representation for the 
in-plane motion, are 

V ni ,n2,n 3 ,m(p) = J dpdzdz' V d {p,z - z') 

(z 1 ). (5) 



x e 



ip-p/h 



For ni = ri2 = na = = we obtain the effective 2D 
interaction ©. Note that although the relative magni- 
tude of a single virtual excitation is as small as gr/ 2 <C 1 
(see App. [Cj, the totality of these processes is essential 
for the correct description of the contribution of short 
distances < Iq (or virtual energies > JvjJq) to the inter- 
particle scattering. We mention also, that within second 
order the p-wave contribution due to the short-range part 
of the interparticle interaction can still be neglected (see 
Ref. [IH for discussion). 



III. APPROACHING THE INTERMEDIATE 
COUPLING REGIME 

The considered problem is characterized by two param- 
eters: g = Pp rd/fr and 77 = Ppho/h. The first parame- 
ter describes the strength of the interparticle interaction 
with respect to the mean kinetic energy: g -C 1 character- 
izes the weakly interacting regime, while we have g > 1 in 



1 This expression is unique up to a momentum independent addi- 
tive constant, which depends on the regularization of the Fourier 
integral at the origin. Such a constant, however, corresponds 
to a short-range interaction and has no physical effect in a sin- 
gle component Fermi gas because its contributions vanish upon 
proper antisymmetrization. Therefore, we set this constant to 
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the regime of intermediate and strong interactions. The 
second parameter 77 describes the strength of the confine- 
ment: 77 <C 1 corresponds to strong confinement. When g 
and r\ are small, g,rj -C f , one can use perturbation the- 
ory to calculate various quantities. Many physical effects, 
however, occur at intermediate values of g (we will always 
assume strong confinement with 77 -C 1), see Fig. [9l for 
which one cannot limit oneself to lowest order diagrams 
or to a specific sequence of diagrams (ladder diagrams in 
a dilute system). To obtain analytic expressions in this 
case, one can use analyticity arguments to extrapolate 
expressions obtained in the weak coupling regime to in- 
termediate coupling strength. Of course, the accuracy 
of such expressions cannot be estimated. However, (pro- 
vided the relevant physics is present in the weak coupling 
regime) they can be used to make qualitative statements 
on the behavior of the system and often obtain reasonable 
quantitative estimates. Following this strategy, one has 
to select a "reasonable" set of Feynman diagrams, which 
allows one to write down a closed set of integral equations 
for relevant physical quantities. This sequence should, of 
course, catch the relevant physics and be consistent with 
general physical principles such as conservation laws and 
particle statistics. For the purposes of this paper we will 
use the (conserving) HFA (see Refs. [U, [32[ and discus- 
sion below), which describes the motion of particles in 
an average potential with exchange effects taken into ac- 
count and, therefore, can be used to describe phenom- 
ena for which interparticle collisions are not important 
(zero sound in our case) . Note that taking exchange con- 
tributions into account is crucial: They guarantee the 
cancellation of all interaction contributions in the case 
of a short-range interparticle interaction (when Vo(p) is 
momentum-independent, Vo(p) = const.), as it should be 
in a single-component Fermi gas. 



IV. QUASI-2D DIPOLAR FERMI LIQUID 



Properties of a normal (non-superfluid) Fermi system 
are conveniently described in terms of Green's functions. 
The single-particle Green's function (see, for example, 
Ref. [1^]) for a two-dimensional system (in the following 
we shall be using units in which K = 1), 



G(p) = -i dt dr e^-P") (T{#(t, r) $ (0, 0)}) , (6) 



where T stands for the time-ordering operator (T- 
product) and p = (w,p), carries information on single- 
particle excitations, while collective behavior resulting 
from two-particle correlations is described by the two- 



particle Green's function 

G 2 ( Pl: p 2 ;p 3 ,p 4 ) = J lJ ( // , ( /r.e" -' •' » r 

i=i 

x (T{ n) $(t 2 , r 2 ) ft fa, r 4 ) ft(t s , r 3 )}) 
= G(pi) G( P2 ) 
x [5{p! -pz)8(]P2 - pa) - S(pi -pi)5{p 2 - P3)] 
+ G(p x ) G( P2 ) r(p u p 2 ;p 3 , P4 ) G( P3 ) G(pa), (7) 

or the closely related vertex function T(pi,p 2 ;p3,p4): 
Single-particle excitations correspond to poles of G(p), 
collective modes and instabilities of the many-body 
system are encoded in poles of the vertex function 
r(pi,p2;P3,P4)- In a homogeneous system, this quan- 
tity depends only on three independent momenta, and 
it is convenient to introduce T(pi,p 2 ;q) = T(pi,p 2 ',pi + 
q,p 2 ~q) such that q is the transferred momentum, which 
satisfies the Bethe-Salpeter equation in the particle-hole 
channel [3!| 

r(pi,P2;g r ) = r P h(pi,P2',q) 

-i { f ph (p 1 ,p' + q:q)G(p' + q)G(p') 



x r(p',p 2 ;g) 



dp' 

(2tt) 3 ' 



(8) 



where r p h denotes the particle-hole irreducible vertex, 
i. e., the sum of connected vertex diagrams with two in- 
coming and two outgoing fermionic lines which cannot 
be divided into two parts by cutting two fermion lines of 
opposite direction. 

The single-particle Green's functions can be expressed 
in terms of the self-energy function through the 

Dyson equation [3^] 

G{p)- l = GW{p)- l -V{p), (9) 

where the non-interacting Green's function is 

GV>{p)- l =u-tb), (10) 

with £(p) = p 2 /2m — [J, and the shifted chemical poten- 
tial fx = fi 1 — u)q/2. S, in turn, is connected with the 
vertex function by the equation of motion (Schwinger- 
Dyson equation) [39| 

E(p) - i J [V Q (p - Pl ) - V (0)} G( Pl )^ 



(11) 



+ / r(p 1 +p2-p,p;p 1 ,P2)G(p 1 )G(p2) 

dpi dp 2 



G(pi +p 2 -p) Vb(p- pa) 



(2tt) 3 (2tt) 3 



The solution of the coupled system of equations ([8]), 
© and pT|) is specified by the irreducible vertex r p h, 
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which is the sum of an infinite set of Feynman diagrams. 
As a result, one cannot write the irreducible vertex T p h in 
a closed form in terms of the Green's function G and the 
vertex T, and some approximation procedure of choos- 
ing a subset of contributions is needed. This procedure 
should be consistent with conservation laws and statis- 
tics of the system. A prescription for generating such 
a conserving approximation is to replace the right-hand 
side (RHS) of Eq. (JTTJ) for £ by a functional of G and 
one-particle propagator is then to be ob- 
tained self-consistently from the approximate equation 
for £ and ([5]), and the approximate irreducible vertex 
T p h which determines T via (|5J), can be found by suit- 
able functional differentiation. In other words, by fix- 
ing an appropriate expression for the self-energy S, one 
uniquely determines the expression for r p h in order to 
make the approximation conserving. 

As discussed in Refs. [^H H2|, the simplest example 
of a conserving approximation taking exchange effects 
into account is the HFA, which we will use in this 
paper. In this approximation, the self-energy is given 
diagrammatically as 



tation 




(12) 



or analytically 



E(p) = i J [V (p p') - Vo(0)] Gtf)-0p, (13) 

and is frequency independent. The corresponding 
particle-hole irreducible vertex is 



pi P2 - q pi P2-q 



pi+q 




pi + q 




(14) 



or 



f ph { Pl ,p 2 ;q) = V (q) - Vb(pi - p 2 + q) 
= f P h(pi P2 + q,q) 



(15) 



and is also frequency independent. Eq. ([5]) then shows 
that T(pi,p2] q) does not depend on the frequencies oj\ 
and uj 2 , i.e., T{p l ,p 2 ;q) = r(pi,p 2 ;g), and we have 

r(pi, p 2 ; q) = f P h(pi P2 + q, q) 

-i f T ph ( Pl - p\<l) G(p' + q)G(p') 



x T(p',p 2 ;g) 



dp' 



(16) 



V. SINGLE-PARTICLE EXCITATIONS IN HFA 

From Eqs. ©, ([TO]), and (US]) we obtain the single- 
particle Green's function in momentum space represen- 



G(w,p) = 



1 



w-e(p) +iOsgn[e(p)]' 
where the quasi-particle dispersion relation is 



e(p) = £(p) + S(p), 



(17) 



(18) 



and the self-energy £, after performing the integration 
over a/ in Eq. (|13p , can be written as [note that the direct 
term does not contribute, as Vb(0) = for the effective 
dipole-dipole interaction 



S(p) 



V (p-p')n(p') 



dp' 

(2^ 



(19) 



In this expression n(p) = 0[— e(p)] is the Fermi-Dirac 
distribution at zero temperature. Equation (|19[) has to be 
solved self-consistently together with the particle number 
equation 



n 2 D 



dp 



n(p). 



(20) 



The problem of finding the solution to Eqs. (119)) and 
is simplified considerably by noting that the evalu- 
ation of the RHS of both equations does not require full 
knowledge of the quasi-particle dispersion relation e(p) 
but only of the Fermi momentum p p , which is determined 
by the requirement that the quasi-particle energy (fT5)) is 
equal to zero for p = pp = pfP, 



£(Pf) = Pp/2m -n + E(pf) = 0. 



(21) 



The solution pp to this equation actually depends on 
the direction p = (cos <fi, sin 4>). However, to shorten the 
notation, in the following we will often write pp and p^ 
instead of pf(4>) and Pf{4>) respectively. 

We proceed by specializing (TK)1) to the Fermi surface 
p = pf and inserting the resulting expression for E(pj?) 
in (EU), 



f^=M+ 1 r dp'p'vV {p F -p'), (22) 

where p' F = pp((f>') and v = m/2w is the density of 
states. An expression for the chemical potential fi can 
be obtained by taking the integral of this formula over 
the angle <f> and making use of the fact that the particle 
number equation (|20[) is equivalent to the condition 



dd> 2 



(0)2 



(23) 



if we express the density as n 2 D = Pp^/^Tt. We find 



_ (o) 1 f dtp dtp 

fX — Ep - 

771 



Eqs. (|22[) and ((24]) form a closed system for the deformed 
Fermi surface pp and we obtain the joint solution to these 
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equations numerically by means of an iterative scheme 
which is described in detail in App. [21 The resulting 
quasi-particle dispersion relation (fT8)l is shown in Fig. [2j 
along with the linear approximation at the Fermi surface 

e(p)nv F (4>)\p-p F (ct>)}. (25) 

Here v F (<p) is the radial component of the Fermi velocity: 

v F (0) =Ve\p F (4>)] 

= pv F ((f)) +e v F (4>). 



(26) 



For future reference we summarize the corresponding 
results obtained in perturbation theory 24]. To first 
order in the dipole-dipole interaction, on the RHS of 
Eq. (fl9l) we insert the distribution function of a non- 



interacting Fermi gas, n(p) = 0(p F — p). The self-energy 
function can then be expressed in terms of complete el- 
liptic integrals. We omit the cumbersome analytical ex- 
pression and content ourselves with stating the result for 

(o) (o) * 
P = P F =Pf P, 

s (D (f (0), 



[f sin 2 0o cos(20) - 2P 2 (cos0 o )] ■ 

(27) 

We write the chemical potential as p = ey + 5p. The 
first order correction 8[i is given by (f2"4"|) [with p' F on the 

RHS replaced by p F ^], 

d >= f34 O) ^(cos0 o ). 



(28) 



Combining these results with (|21l) we find the equilibrium 
deformation of the Fermi surface 5p F {4>) = p F {4>) — p F , 



p F 



(29) 



^9P F sm 2 9 Q cos(2^>). 



From its definition in (|26[) , the radial component of the 
Fermi velocity is then 



A o 

10 



^ 2 



10 



0.5 



1.5 



P/P 



CO) 




0.5 1 1.5 2 

p/pj? 

FIG. 2. (Color online) Quasi-particle dispersion relation for 
g = 1, 9o = tt/4 (top) and g = 2, O = tt/2 (bottom) in HFA 
(solid lines) and first order perturbation theory (dashed lines) 
for some values of <j> p [p = p (cos <j> v , sin 4>p)\ (lowermost blue 
curves: 4> p — 0, middle green curves: 4> p — tt/4, topmost 
red curves: <f}p — 7r /2). The linear approximation (|25p at the 
Fermi surface is shown as dotted lines. Note that the results 
obtained in HFA and perturbation theory agree very well for 
g = 1, 8q — 7r/4, whereas deviations are clearly visible for 
9 = 2, O = 71-/2. 



v F (4>) 



,(°) 



[1 + ^ ffj p 2 (cos(? ) _ ^ 5 sin 2 # 008(20)] , 

(30) 



Finally, the (radial) effective mass is defined as m*(<p) 
p F ((f>)/v F (<f>). Then, for the deviation 8m((f>) — m*((f>) 
m of the effective mass from the bare mass we have 



5m((f>) 



d£«(p(?) 



pf 



dp 



(31) 



= -±gP 2 (cos6 ) + ^«?sin 2 o cos(20). 



VI. 



COLLECTIVE EXCITATIONS IN HFA: 
ZERO SOUND 



Collective modes u> = w(q) of the system correspond 
to poles of r(pi,p2;q) with respect to the variable w, 



which is the frequency component of q. In the vicinity of 
a pole we have T 3> r p h and, therefore, we may neglect 
the first term on the RHS of Eq. (TIB)) . In the result- 
ing homogeneous equation the second argument p 2 of T 
acts as a parameter. Hence, near its pole, the function 
T can be represented as a product x(Pli l) x'(P2! q) 01 
two functions. After cancelling x'(P2! q) on both sides of 
Eq. (|16[) and integrating over lu' , we obtain 



X(p';q). (32) 



f dp' ~ , . 

x(p; q) = J r ph (p p , q) 



(2tt)' 
n(p') - n(p' + q) 



cj + e(p') - e(p' + q) + zO sgn(w) 



The quantity e(p' + q) — e(p') on the RHS of the last 
equation is just the energy cost of creating a particle- 
hole pair by exciting a particle from a state p' within 
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the Fermi surface to a state p' + q outside the Fermi 
surface. Therefore, a stable collective mode is possible 
only when the energy (frequency) uj of the mode lies out- 
side the particle-hole continuum (PHC) (so we can omit 
the imaginary term in the denominator). Otherwise the 
integrand has a pole at u> = e(p' + q) — s(p'), which ulti- 
mately leads to strong Landau damping of the collective 
mode [39(. 

In the long wavelength limit |q| — > 0, the main con- 
tribution to the integral (j32|) comes from states in the 
vicinity of the Fermi surface, and we can rewrite this 
equation as 

, s f #' f, / t ns p' F l m q ' V F i i s 

(33) 

where x(p) = x(p;0), p' F = Pf (<£')> and v F = vp (<£')• 
The function x(p) is thus completely determined by its 
values on the Fermi surface and by setting p = pf we 
can obtain an equation for the restriction of x(p) to the 
Fermi surface, 



x(Pf) 



where we define the dimensionless quasi-particle interac- 
tion function (/-function) as 



i^,0')^r ph ( PF -p^,o) 

= v[V {0)- Vo(pf-Pf)] 
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|P_F ~ P F \ 

(0) 

Pp 

P 2 (cos9 )- ^u(p F - p' F )sin 2 6 



(35) 



Note that only the exchange interaction contributes to 
the /-function. For 6q = 0, i.e., when the dipoles are 
perpendicular to the ccj/-plane and the system is symmet- 
ric with respect to rotations around the z-axis, the Fermi 

momentum is isotropic and equal to p F ^ [see Eq. (|23[) ]. 
and (|35p simplifies to 



F(cf>,(t>') = 2gsm(\<f>-<p'\/2), 9 



(36) 



Equation Q34p shows that uj depends linearly on q = 



(o) 

uj = v F s q. 



(37) 



Due to the anisotropy of the dipole-dipole interaction, s 
will in general be a function of the propagation direc- 
tion (j)q [q = q (cos<^ q ,sin0q)]. The symmetry of the 
problem, however, requires the dependence of s on q to 
be 7r-periodic and even, hence, it is sufficient to restrict 
ourselves to the range < <^ q < n/2. 

As it has already been pointed out, the excitation en- 
ergy oj of the collective mode has to be separated from the 



PHC. Equation (|3"4"]l shows that in the long wavelength 
limit q — > this requirement reduces to 



,(°) 



s > u ph , 



(38) 



where the quantity u p h is the slope of the upper boundary 
of the PHC in the direction q at q = 0, and can be 
computed by taking the maximum over all values of the 
angle <f)', 



Vph = max{q ■ v F }. 



(39) 



For the purpose of solving Eq. (j34]) numerically it is 
convenient to replace x by another function 



v(4>) 



q • v F 



P V F ^ 0) s -q. VF 



x(pf)- 



(40) 



Equation (|3"4")l then becomes 



(v^'s-q- v F )v(4>) = 



1 q-VF 
m p • vf 



2tt 



(41) 

One should note that Eq. (|4"Tj) - in contrast to Eq. (jM)) 
- allows for non-trivial solutions even for g = [or 
F(4>,4>') = 0]. These solutions with s e [-1, 1] take the 
form of (5-function in <p and correspond to single particle- 
hole pairs from the PHC. 

We solve Eq. (|4ip numerically by discretizing it on an 
evenly spaced grid of 2000 points in the variable </> and 
rewriting the integral on the RHS as a sum according to 
the trapezoidal integration rule. Convergence tests show 
that increasing the number of grid points further by a 
factor of 2 results in a change in s that is less than 10~ 5 
(for all values of the parameters in the range 0.1 < g < 2, 
and < #o, 0q < 7r/2; going to g <C 1 is not possible in 
this approach, see discussion below) . 

Our results for g — 1 are shown in Figs. [21 IH and [S] 
We see that the existence of zero sound and the value of 
the sound velocity strongly depend on the propagation 
direction, on the tilting angle, and on the strength of the 
interaction. There is no dissipationless zero sound mode 
if the tilting angle is smaller than some critical value, 
see Fig. [31 as then all numerically calculated eigenval- 
ues s are below or equal to the limiting velocity w p h of 
the PHC, thus violating the propagation criterion f|38|) . 
This is to some extend similar to the propagation of zero 
sound in the homogeneous three dimensional case, which 
was studied in Refs. [24], |33[: The authors of these ref- 
erences found that there is no undamped propagation of 
zero sound for a wide range of angles perpendicular to 
the direction of the polarization of dipoles. The reason 
for this at first glance counterintuitive statement is that 
the contribution of the direct interaction (which is repul- 
sive in coordinate space) to f p h [or to F(cf), 4>')] vanishes, 
see Eqs. (IT5|) and (@}, such that the long wavelength col- 
lective behavior of the dipolar gas is governed by the 
exchange interaction (the /-function contains only the 
exchange contribution). (The collective modes without 
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0.025 



0.015 




tt/4 37r/8 tt/2 



FIG. 3. The condition (|38|) for the propagation of zero sound 
is met in regions I and II, where the gray levels encode the 
difference Vp's — v p h for the largest eigenvalue s of Eq. (|4Tjl . 



the exchange contribution were considered in Ref. [29j, 
and, as a consequence, their result is determined by the 
momentum-independent term in the Fourier transform 
Vo(q) of the dipole-dipole interaction which is omitted in 
our paper because it is cancelled by the corresponding 
exchange contribution.) 

We first find a zero sound mode that satisfies the prop- 
agation criterion ((38|) (see Figs. [3] and |H top) for 6*o ~ 0.5 
and q = (i. e., in the direction of the projection of the 
dipoles d on the xy-plane). It is a longitudinal mode that 
is concentrated symmetrically around its propagation di- 
rection. As 6q increases, this mode can propagate in a 
broader range of angles up to </> q ps 7t/4 (region I in Fig. [3] 
see also Fig. O top); However, the corresponding sound 
velocity drops below the propagation boundary Vph at 
O « 0.9. 

A different mode emerges from the continuum at 9q ss 
7r/4 and <7J q = ir/2 (region II in Fig. This mode is 
antisymmetrically peaked around the direction of prop- 
agation. At even higher values of 9q we find more than 
one modes that satisfy (f38|) (see Figs. HI bottom, and [5j 
bottom). 

This behavior of collective modes remains the same 
when we increase the interaction strength to g ~ 2. In 
particular, the angles 0q w 0.5 and 9q s» tt/4 at which 
the symmetric and antisymmetric modes appear, respec- 
tively, are practically left unchanged. However, we find 
that the zero sound modes become more "distinct", i.e., 
the curves in Figs. 0] and [S] are separated further from the 
PHC. 

For g < 1 the above-mentioned peaking of modes 
around the forward direction is even more pronounced, 
and a high number of grid points is required to prop- 
erly resolve these modes, making it impossible to go to 
very small values g <C 1 with our numerical method. 
At 0.1 < g < 1 it becomes increasingly difficult to 
make quantitative statements as with regard to the re- 




tt/8 



tt/4 



3tt/8 



7r/2 



FIG. 4. (Color online) The speed of zero sound as a function 
of the tilting angle 8o for g = 1 and q = 0, tt/2 (top, bottom). 
Note that for (frq = tt/2 and 8o close to tt/2 we have two zero 
sound modes. The region that is forbidden by the condition 
((38)) is shaded. 



gion where the condition (|38p met, since the quantity of 
interest Vp^s — v p h is of the same order of magnitude as 
its estimated error. Qualitatively, however, we find the 
same behavior as described above for 1 < g. 

We also solve Eq. (|52")l with finite values of q. Some 
results are shown in Fig. [5] and demonstrate that the 
dispersion uj = w(q) of the collective mode obtained in 
this manner agrees well with the linear approximation 
from Eq. (|37|) . For details of the numerical procedure see 
A PP . El 



VII. INSTABILITY OF A SPATIALLY 
HOMOGENEOUS SYSTEM 

Apart from poles on the real axis corresponding to 
excitation frequencies of collective modes, Eq. (1521 can 
also have purely imaginary eigenvalues that indicate the 
existence of unstable modes growing exponentially with 
time. At the onset of an instability we have an eigenvalue 
uj = 0, and, depending on the corresponding value of the 
wave vector q = |q|, one has a long wavelength instabil- 
ity (q — 0) resulting in local collapse of the system, or 
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tt/8 tt/4 3tt/8 tt/2 



0.0 



0.8 




FIG. 5. (Color online) The speed of zero sound as a function 
of the propagation direction <^> q for g — 1 and #o = 7r/4, 7r/2 
(top, bottom). Note that for 9o = 7r/2 and </> q > 7r/4 we 
have two zero sound modes. As in Fig. [3] the region that is 
forbidden by the condition (|38[) is shaded. 




FIG. 6. (Color online) Collective mode dispersion relation for 
<j!> q = 7r/2, g = 2 and So = 7r/2 (solid lines). For comparison 
we also plot the linear approximation (I37|) as dashed lines. 
The inset shows that at long wavelengths we have two zero 
sound modes. 



finite wavelength instability (q ^ 0) leading to breaking 
of translational invariance and formation of a periodic 
spatial structure (density waves). 



A. Long wavelength instability 

Equation (j4"Tj) shows that, if the instability occurs at 
q — > 0, the boundary of the instability region is deter- 
mined by 



l 



mvF(4>) J 27r 



F(q>,fi)p F (cj>')vW), (42) 



and, hence, the instabilities could occur only in the 
regime of intermediate coupling, g w 10 

Equation (|42l) is equivalent to the Pomeranchuk cri- 
terion [34] on the Landau /-function to ensure stability 
of a three-dimensional isotropic Fermi liquid [39]. We 
briefly review this method and its generalization to the 
two-dimensional anisotropic case. 

In the framework of Fermi liquid theory [401 ] . the 
change of the quasi-particle momentum distribution 
function 5n(p) results in the change in the energy density 



SS 



dp 

! L 

1 f dp dp' 



(2tt) 2 ~ ^ 5ni ^ 



(2tt) 2 (2tt) 2 



f(p,p')Sn(p)6n(p'), (43) 



where f(p,p') is the quasi-particle interaction function 
(/-function). A slight distortion Spf(4>) [note that in 
this section, 5pf is used another way than in Sees. fVl and 
IVIII] of the anisotropic Fermi surface pf(4>) corresponds 
to a distribution function of the form 

Sn(p) = e\p F {4>)+SpF{<f>)-p] - 0\pf{4>) -P\. (44) 

Following Pomeranchuk [34j we expand 8£ in Spf'- The 
first order term vanishes due to the fact that the energy 
of the Fermi liquid, considered as a functional of Spf, is 
stationary at 8pp = 0. To second order we have 



68 



Attiti 2 



2tt 



p F {(j>)m v F (</>) Sp F (</>) 2 



+ / ^^PF{4>)PF^)F^A')5pF^)5p F ^') 

(45) 

where F(<fi,<f)') = v /(pfjPf) i s ^ e dimensionless quasi- 
particle interaction function already introduced in ([35]) • 



2 In Eq. 14H , on the other hand, the prefactor on the LHS can 
become arbitrarily small when Vp s gets close to the limiting 
velocity of the PHC v p ^, and there is no generic restriction on 
the size of g for this equation to have a solution, i.e., for zero 
sound to occur. 
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Thermodynamic stability requires that the variation of 
the energy density be positive for any choice of Sp F . This 
condition can conveniently be formulated by expand- 
ing the 0-dependent functions in Fourier (double) series 
f{4>) = E m fm e™ , where f{4>) is one of p F {<j>), 8p F {4>), 
and v F (4>), and F(4>,<j/) = ^ Fm,m> e*™*-™' +"> , 
where the coefficients F m .m' take the role of generalized 
Landau parameters. 

Replacing the functions in the expression for 8£ by 
their respective Fourier series expansions we obtain 

S£ = ^ Sp* m M m . m , 5p m >, (46) 

m,m' 

where the entries of the (self- adjoint) matrix M are 

Pm—k v m'—k + 'J2Pm-kFk,lPm'-l- ( 47 ) 

A 2D anisotropic Fermi liquid, therefore, is thermody- 
namically stable if and only if all eigenvalues of M are 
positive. In the case of isotropic interactions, we have 
M mjm > — 5 m ,m' Mm, and the requirement of stability re- 
duces to the usual Pomeranchuk criteria, i.e., M m > 
for all to € Z. 

At the onset of an instability there exists an eigenvalue 
that is equal to zero, i.e., there is a non-trivial solution 
to the equation 

M ™,m> 5pm- = 0. (48) 

m' 

In terms of the original Fermi surface deformation Sp F 
this equation reads 

Sp F (4>) + ^T-T / ^-F(cf>, ( f>')p F ( ( f>')Sp F (cf>')=0, 

and if we identify Sp F with v we are lead back to Eq. (1421) . 
The line in the g6>o-plane on which this equation has a 
non-trivial solution determines the boundary to the red 
shaded region that is labeled LWI in Fig. |H1 In the in- 
terior of this region there exists an unstable long wave- 
length mode. 

B. Finite wavelength (density wave) instability 

An instability at finite momentum q = |q| drives the 
system towards a state with stationary fluctuations in 
the density. We address this problem numerically (see 
details in Appendix B). After solving Eq. (|B3|) [which is 
equivalent to Eq. Q32p] numerically for q ^ 0, we find 
such an instability with </> q = tt/2 (i. e., q is perpendicu- 
lar to the x-axis - the projection of the dipoles d on the 
xy-plane) and q rs 2p F (<j> tl ) in the blue shaded region 
in Fig. IH1 As can be seen from Eq. (|B3[) , the eigenvec- 
tors ^q(p) oc (ap +q / 2 a p-q/2) which correspond to an 
instability, signal the formation of a density modulation 



with momentum q (density wave). On the other hand, 
the dependence of ^ q (p) on p excludes a description of 
the transition in terms of a simple local order parame- 
ter (r) ip(r)) = p(q)cos(q- r). In the isotropic case 
with #o = 0, where there is no preferred direction and 
the system is invariant with respect to rotations around 
the z-axis, the instability occurs at g sa 1.45 and is in- 
dependent of the angle q . This result is in agreement 
with the value g ~ 1.42 which was found by the authors 
of [13] ■ Previous studies within the RPA [25], [3(| pre- 
dicted a considerably smaller value of g = 0.5, which is 
due to the fact that the RPA overestimates the effects 
of the interparicle interaction as it neglects the exchange 
contribution. 

Note that for r p h(p — p',q) = Vb(q) (i-e., when only 
the direct interaction is taken into account and the ex- 
change one is neglected), x(Pi Q) i s p- independent and 
Eq. ([32D reduces to 

l-yo(q)n<°>(w,q) = 0, (50) 

where 

p' + q 

n(°)( w , q )=<^> 

p' 

= f dp' n(p')-n(p' + q) 

J (2tt) 2 u + e(p') - e(p' + q) + iO sgn(cj) 

(51) 

is the 2D polarization operator. Equation (|50p is used 
to study long wavelength (q — > 0) plasmon oscillations 
in electrically charged systems (see, for example, [41]). 
Although keeping only the direct interaction in the long 
wave-length limit is legitimate for Coulomb systems (be- 
cause of the divergence of the Coulomb interaction for 
a small transferred momentum q — > 0, while the ex- 
change one is finite due to a non-zero momentum trans- 
fer, |p — p'| ~ Pf), this approximation gives physically 
incorrect results in a Fermi system with a finite Fourier 
transform of the interparticle interaction for small mo- 
mentum transfer (like in the considered case of a dipolar 
monolayer). In this case, the direct and the exchange 
contributions are of the same order and keeping only the 
former yields unphysical results. For a short range in- 
terparticle interaction (with a momentum-independent 
Fourier transform), the two contributions have to can- 
cel each other, resulting in no interaction effects in a 
single-component Fermi gas with a short-range interac- 
tion. Similar considerations also apply to the analysis of 
instabilities in a dipolar system on the basis of Eq. (|32p : 
Keeping the exchange contribution in this equation is es- 
sential in order to obtain correct results consistent with 
fermionic statistics of particles. 

In order to take the exchange contribution into account 
in Eq. (|50[) . one can modify the polarization propagator 
by including an entire interaction ladder in the polariza- 
tion bubble (this is equivalent to including the exchange 
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interaction in f p h), see also Ref. [37 



n( w ,q) = n(°)( w , q )+<^]3> 



■<xx>- 



(52) 



(Another possibility is to include a local field correc- 
tion [35l 1 3 61] , similar to the consideration of the density- 
density response in Coulomb systems. This approach, 
however, is very sensitive to a particular choice of the lo- 
cal field correction.) With the modified polarization op- 
erator n(u;,q), the density-density correlation function 
(in frequency-momentum space) reads 



X(w,q) 



n(w, q ) 



i-y (q)n(w,q)- 



(53) 



and Eq. (|50l) with n(°)(u;,q) replaced with n(w,q) cor- 
responds to the instability in the density-density corre- 
lation function. Note that when Vb(q) is replaced by a 
momentum independent constant V$ (that corresponds to 
a short-range interaction), the modified polarization op- 
erator is n(w,q) = n(°)(w,q) [l + y n(°)(w,q)] _1 , and 
the density-density correlation function reduces to the 
polarization operator of a non-interacting gas, %(w,q) = 
n(°)(o;, q), as it should be in a single-component Fermi 
gas. 



VIII. SUPERFLUID TRANSITION 

We now discuss the superfluid instability in a dipolar 
monolayer. As we will show, this instability is sensitive 
to the details of two particle scattering and, as it was 
already pointed out at the end of Sec. HH this requires 
us to take into account contributions of short distances 
< lo and high energies > ojq. Therefore, in the follow- 
ing we will not limit ourselves a priori to configurations 
with all particles residing in the ground state of the trap- 
ping potential but rather allow for virtual transitions to 
arbitrarily highly excited states. 



A. Gap equation 

The superfluid transition is characterized by the order 
parameter (gap) A(r,r') = Vd(r — r')(i/>(r) ip(r')), which 
attains non-zero values for temperatures below the crit- 
ical temperature T c . Note that in two dimensions at fi- 
nite temperatures long-range order is actually destroyed 
by phase fluctuations and the mean-field order param- 
eter is zero. The superfluid density, however, remains 
finite and the transition to the superfluid phase follows 
the Berezinskii-Kosterlitz-Thouless scenario |42l - |44j and 
occurs at a temperature y c ( BKT \ Nevertheless, we may 
consider the mean-field critical temperature T c as a re- 
liable estimate of the value of y c ( BKT ) as the difference 




between the two of them is small in the weak coupling 
regime (45j . 

For < T c — T -c T c , the order parameter is a solu- 
tion to the homogeneous Bethe-Salpeter equation in the 
Cooper channel or linearized gap equation [39j . 



(54) 



where r pp denotes the particle-particle irreducible ver- 
tex, which is the sum of connected diagrams with two in- 
coming and two outgoing lines which cannot be divided 
into two parts by cutting two fermion lines of the same 
direction. 

In these diagrams, thick lines correspond to interacting 
Matsubara Green's functions in the presence of the trans- 
verse trapping potential. Following Gor'kov and Melik- 
Barkhudarov [38[ we expand the RHS of Eq. ((54)) to sec- 
ond order in the dipole-dipole interaction Vd, 




Here, with each thin line is associated a non-interacting 
Matsubara Green's function which can be written in the 
form [H 



G(oj s ,p,z,z') = ^ c t ) n{z)cj) n {z')Q n (uj S: -p) 1 



n=0 



(56) 



g n (uj s ,p) 



iu s - £(p) - oj n 



where £(p) = p 2 /2m — /1, and fermionic Matsubara fre- 
quencies are u) s = (2s + 1)ttT for integer s. The chemical 
potential is fj, = Ep -\-Sfx, where the first order correction 
8[i is given by Eq. ([28| [due to the exponential smallness 
of T c (see discussion below) we may use the zero temper- 
ature value of //]. Thus, strictly speaking (1551 contains 
terms of infinite order. It is, however, convenient to per- 
form the expansion in S/i at a later stage. 



Of the diagrams on the RHS of (|55|) the first one gives 
the leading (first order) contribution, the others corre- 
spond to (second order) corrections: The second diagram 
is obtained from the first one by inserting an exchange 
self-energy part (note that the direct term is absent for 
the dipole-dipole interaction, see Sec. [Vj). This diagram 
comes with a factor of 2, as we could have equally well 
inserted the self-energy part in the upper particle line. 
The quantity SVd represents second order corrections to 
the bare dipole-dipole interaction and is given by the set 
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of diagrams 



sv d = sv} a) 



sv^ = 



sv 



(*>) 



sv 



sv 



(d) 



(57) 



sv 



(c) _ 





sv^ = 



These diagrams describe processes in which one of the 
incoming particles polarizes the medium by exciting a 
virtual particle-hole pair. In diagram (a) the particle 
and hole annihilate each other while interacting with the 
other incoming particle, whereas in (b, c, d) the hole is an- 
nihilated by one of the incoming particles: by the second 

I 



incoming particle in (6, c), and by the very same particle 
that created the particle-hole pair in the first place in 

In writing down the analytical expressions that cor- 
respond to the diagrams in (|55[) it is convenient to la- 
bel particle lines by two-dimensional momentum vectors 
P = {pxiPy) and HO quantum numbers n, since the Mat- 
subara Green's function f|56[) is diagonal in this basis. 
Accordingly, we expand the order parameter as 



A{p, z,z') = 



<\>n{z) dJn'iz') A n , n '(p). 



(58) 



n,n'— 



The matrix elements of the dipole-dipole interaction are 
defined in Eq. and we employ an analogous definition 
for matrix elements SV ni ,n 2 ,n 3 ,n4, of SVd- This allows us 
to write the diagrams on the RHS of (f5"5"|) as 



J£) = - I 53 ^«i 1 »2,™3,"4(p-pO T ^^(^'P')^(- w SJ -p') A n3,n 4 (p')7^ 

\ Sj£) = ~ 2 53 Kn,n 2 ,n s ,ru(P-p') 71 X]^3(' Js ' p/ )^ n *( _a;< " _p/ ) 

( 713,714,715 s 

T 53 Kn,n 2 ,n 3 ,n 3 (p - p') Gn 3 (U) a , p') , 
53 ^"2 ,n 3 ,n 4 (P - P') T 53 ^"3 > P') £»4 (~W S , -p') A n3i „ 4 (P')?2^2 • 



m dp' 




(59) 



(60) 
(61) 

(62) 



713,714 



The sum over Matsubara frequencies in (|59p and ([52 
can be evaluated by rewriting it as a contour integral [4f 
and gives 

£n,n'(p) = T 53^ n ( ti;; ''' P ) ^n'(~ w s,-p) 



tanh 



[ j(p)+"n" 
L 2T 



tanh 



C(p)+"o»' 
2T 



2[2£(p)+w (n + n')] 



(63) 



For n, n' ^ and in the quasi-2D limit » /i the 
denominator in (|63[) is positive for all values of p and 
the arguments of the hyperbolic tangent functions are as 
large as ujq/T. Then, to within exponential accuracy in 
this ratio, we may approximate 



£n,ra'(p) 



1 



2£(p)+w (n + n')' 



(64) 



If only one of n,n' is different from zero, the denomi- 
nator in (jSB"]) is still positive for all p. For concreteness, 



I 

let us assume that n ^ and n' = 0. Then we set 

1 + sgn[£(p)] 



£n,o(p) 



2[2£(p)+w n] 



(65) 



thereby omitting a subdominant contribution to the in- 
tegral (|59p from a narrow shell in momentum space 
l£(p')l ^ ^ where tanh[£(p)/2T] deviates appreciably 
from sgn[£(p)]. 

Finally, in the case n = n' = 0, (|6"3")l reduces to 



£[£(p)] = /C ,o(p) = tanhK(p)/2T]/2£(p) 



(66) 



In this expression, the denominator vanishes on the Fermi 
surface, i. e., for £(p) = 0. This leads in the limit T//x <C 
1 to a logarithmically divergent integral cx ln(T//x) on 
the RHS of U5HD when n 3 = n 4 = 0. 

We have to keep contributions of second order in the 
gap equation only if they are multiplied by such a large 
logarithm |38j |: In (1601) this is the term with 713 = 714 = 
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ri5 = which contains the factor 



IdK, 



T p) Go(-u> a , -p) 2 = - ^K(p)], (67) 



and in f|62|) it is the summand with 713 = ri4 =0. These 
simplifications allow us to rewrite Eqs. ([55]) . (|S(J)) , and 

© as [A(p) = A ,o(p), S (1) (P) = 4ld(p)] 



" ^ 41PJ (2tt) 2 



Jt) = - ^2 V nun2 , n3 , ni {p-p')}C n3 , ni (p')A 
j^B = - / IWoCP " p') f K(p')] ^ (1) (P') A( P ')^ 




<5K 



,o,o(p-p')^K(p')]A(p 



d P' 



(68) 
(69) 
(70) 



where JC n , n '(p) and /C[£(p)] are given by Eqs. (l64|) and 
(|55|) , and dni]) respectively. 



B. Renormalization 

Apart from the region p' ss p§\ major contributions to 
the integral in (|68|) that are actually divergent and need 
to be cut off come also from high momenta. This re- 
gion, however, is related to short interparticle distances, 
at which the presence of other particles becomes irrele- 
vant and the dynamics corresponds to the scattering of 
just two particles in vacuum. In the gap equation this 
two-body physics can be taken into account by expressing 
the bare dipole-dipole interaction in terms of the vertex 
function 7 ni ,n 2 ,ri3,n 4 P, p') for two particles in vacuum 
with a total energy E in their center-of-mass reference 
frame. The Lippmann-Schwinger equation for the vertex 
function |47j can be represented diagrammatically as 





(71) 



where particle lines correspond to non-interacting zero 
temperature Green's functions, 



g<°>( Wj P , z, ■/) = J2 M*) M*') G i 0) (^ p). 

n—0 

Gf( Wj p) 



(72) 



u> — p 2 /2m — ljqu + iO 



r 



For our purposes it is convenient to rearrange the order 
of terms on the RHS of (I7T1) as 



7m, 



ri3,"4 C^i P' P ) — ^ni ,ri2,ri3,ri4 (P P ) 

/ r yn 1 ,n 2 .n 5 .n s {E,p,CI L )K n ^ n6 {E,q L ) 

n s ,n 6 

dq 



X Vn s ,r. 



L(q-p') 



(73) 



where the kernel K n ^ n r(E, p) is given by the integral over 
frequencies, 



K n>n ,(E,p) = -i J G^(E + io,p)G ( t ° ) (-u J ,-p) 



p 2 /m + u)o(n + n') — E — iO 



(lid 

2^ 



(74) 



Below we shall choose E = 2/i, in which case ([741) be- 
comes 

1 



2f(p) +w (n + ri) - iO' 



(75) 



In order to carry out the renormalization in the most 
transparent way we rewrite the gap equation (I55I) and 
the Lippmann-Schwinger equation (I73I) schematically as 

A = -VKA-VdKT,^A-8VKA, (76) 
7 = (1 - jK)V. (77) 



We "multiply" ([76]) from the left by 1 - jK. Then, using 
(1771) and neglecting terms that contain 7VEW and ^/SV 
(those are contributions of third order), we obtain the 
renormalized gap equation 



"7 



(AC - K)A - VBJCEWA - 5VKA. 



(78) 
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In the renormalized gap equation, the kernel /C n , n ' (p) 
has been replaced by the difference /C„, n '(p) — 
K n<n i (2/i, p) which converges rapidly at high energies, 
i.e., high p and large HO quantum numbers n and n' . 
To wit, for n = n' = 0, from Eq. (|66p we have 



£o,o(p) - #o,o(2/i,p) = 



tanh[£(p)/2T] 
2?(P) 



1 



2e(p)-*o' 

(79) 

which decays as 1/p for p — >• oo, ensuring the conver- 
gence of the integral over the momentum in the first term 
on the RHS of (|78l) without the need to introduce an ad- 
ditional cut-off. If either n or n' is nonzero, with (|65l) we 
have 



£n,o(p) - K nfi (2fj,,p) PS - 



K(p) 



2£(p) + uj a n' 



(80) 



In l|78p this results in a term that is 0(/i/wo) and may 
safely be neglected. Finally, for both n and n' not equal 
to zero, with the aid of (IM|) we obtain 



£n,n> (P) - Kn,n> (2/i) p) ~ 0. 



(81) 



Therefore, in the propagators in (|78[) . we may re- 
strict ourselves to the HO quantum numbers being 
equal to zero, and we obtain a closed equation for 
A(p) = A 0:0 (p): With j(E,p,p') = 70,0,0,0 (-B,P,P') 
and 8V (p, p') = £>Vb, ,o,o(p, p') we have 



A( P ) = - y| 7 (2M,P,p') 



jc(0 - 



i 



2f - *0 



9/C 



vb(p- P ')s (1) (p')^(e') 
+ «Vo(p,p')^(o1a(p' 



dp 7 



(82) 



To proceed we need to find an expression for the vertex 
function 7(2/t, p, p'). 



C. Two-body vertex function 



Iteration of Eq. (|73|) yields the familiar Born series. 
Terminating this series at second order we obtain for 

n\ = n-2 = 123 = = 

7 (2/i, p, p') = 7 (1) (2/i, P, P') + 7 (2) (2/*, P, p'), (83) 

where the first order contribution is just the Fourier 
transform of the effective 2D dipole-dipole interaction 



7W(2/i,p ) p , )=Vo(p-p / ), 



(84) 



and to second order we have to include virtual excita- 
tions to higher HO levels at the intermediate stage of the 
interaction, i. e., we have to take the sum over n, n' £ No, 



7 (2) (2/ / ,p,p') 

-E 



Vn,n< (P - q) Vn,n< (q ~ P') dq 

2£(q)+wo(n + n')-iO (2tt) 2 



(85) 



where V n , n '(p) = Vo,o,n,n'(p) = V n ,n>,o,o(p)- Explicit 
expressions for these matrix elements - calculated exactly 
as well as in WKB approximation - are given in App. [C] 



D. Asymptotic gap equation 

In the integral on the RHS of Eq. (1521 we perform the 
change of variables 



P = p(£» = VMC 



4 0) )p, 



(86) 



where p = (cos (j>, sin <f>). We shall simplify the notation 
by introducing the functions (avoiding to explicitly state 
the dependence on £ and <j)) 

f{0 = J ^[7(2M,P,P , ) + W>(P,P')] A(£',^) 

g(e) = J ^u^(2n,p,p f )A(e,n 
na = J f^ (p-pos (i) (p')A(r,n 



(87) 



With these definitions we can write the gap equation 
in a very compact form, 



A = - 



rf£' 



£(0/(0 + K'(OMO- 



In the second term on the RHS we integrate by parts. 
Neglecting boundary terms that are of second order in 
the dipole-dipole interaction and are not multiplied by 
the large logarithm ln(T//t), we obtain 



A = 



n I 2£ - tO 



(89) 

In the limit T/fi -C 1 the function /C(£') behaves as 
/C(£') ~ 1/2|£'|. Hence, the main contribution to the 
integral over the first term on the RHS of Eq. (|8"9"|) . as has 
already been discussed above, is logarithmic in T j \i and 
comes from states near the Fermi surface where |£'| <C /i. 
In order to single out this contribution we divide the 
integral over £' into two parts: (a) the integration of 

[/(°) - h '(°)} from "A 1 t0 and ( b ) the sum 01 
the integrals of £(£') [/(£') - h'(g) - /(0) + fc'(0)] over 
[-/i, M ) and of £(£') [/(£') - over (/t, oo). In part 

(a) we use the asymptotic formula 



" , tanh(£72r) / 2e> 
-ju 2£' ~ ^ 7tT 



(90) 



where 7 ss 0.5772 is the Euler constant, and in part (b) 
we once more [cf. Eq. (1651) ] replace tanh(£'/2T) by sgn(£') 
and integrate by parts. Consequently, keeping terms that 
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are 0(g 2 ) only when they are multiplied by the large The value £ = corresponds to the momentum 
logarithm ln(T///), Eq. ([551 takes the form, 



-In 

1 

2 



2eT/A 



[/(o)-fc'(o)] 



In 



df' In 



%[/(e')-^'(0]+^5(o). (9i) 



= = ^ 2m(ef + 5/i) « p^ 0) + -^U/z, (92) 



(o) 



where <5/x is the first order correction to the chemical po- 
tential (|28|). Inserting this expansion as well as the ex- 
plicit expressions (|87p for /, g and /i in the gap equation 
JTT) we obtain 



A(£, 



l n | £L§L 



2tt 



^o(p-pf) 



m ay (p- P y o)/ ) 

(o) 
P F 



^ + 7 (2) (24 0) ,p,p^) 



(o)\ 



9p* 



+ ^ (p,p^') 



A(0, 0') - — [Vb(p - p') £« (p') A(£', 0') 



In 



.(0) 



J ^^ o (p- p ')A(£',0')] ^^o(p-P^ 0) ')A(0>')- (93) 



Our calculation of the critical temperature is carried 
out in two steps: First we omit all terms on the RHS 
of Eq. but the first one. This is equivalent to the 
commonly used BCS approach and allows us to obtain 
the controlling factor for the dependence of T c on g as 
well as the leading behavior of A(£, 4>) in the limit g <C 1. 
We will find that T c oc Bp exp(l/Ao), where Ao oc g. The 
calculation of the pref actor of exp(l/Ao) in the second 
step requires us to take account of all terms on the RHS 
of Eq. , which is referred to as the GM approach. 



E. Critical temperature in the BCS approach 

Keeping only the dominant contribution on the RHS 
of Eq. (|9"5]l. we have 

A(^)=ln(-J^ J ^^ o (p-p^')A(0>')- (94) 

We expand A(0, 4>) in terms of a complete set of eigen- 
functions $ s (0), s € No of the integral operator with 

the kernel Vb(p^? — p^ ), which for Ppho -C 1 can be 
written as [cf. Eq. (gj)] 

Voip^ - pf) « -{9/v) sin(|^ - <t>'\/2) 

x [cos(0 + (j)') sin 2 6» + 2 P 2 (cos 9 )] . (95) 

The expansion of A(0, (p) reads 



A(O» = ^A S $ S (0), 



(96) 



where the functions $ s (0) satisfy the eigenvalue equation 
/ ^vV { V f -vf')<5> s {4>') = \ s <5> s {<p), (97) 

and arc normalized to unity according to 
f(d<j>/2n) ^ s ((j)) 2 = 1. We label the eigenvalues A s 
such that A s < A S ' for s < s'. Inserting the expansion 
([9"6")l in Eq. and specifying the resulting equation to 
£ = we obtain the set of equations, for s E No, 



1 - A s ln(T/e£ 



(oh 



0. 



(98) 



Thus the existence of a non-trivial solution for A(0, <p) 
requires at least the smallest eigenvalue Ao to be nega- 
tive [note that ln(T/e^ 0) ) < since T/e^ < 1]. Then 
the controlling factor of the critical temperature follows 
immediately from 1 — Ao ln(T c /e^?') = 0, and we have 
Ao 7^ whereas A s = for s ^ 0. 

Solving the eigenvalue equation (|97p numerically re- 
veals that, as the tilting angle 0q of the dipoles is in- 
creased, Ao becomes negative at 8 C w 0.72, and for 8q not 
too close to 8 C the corresponding eigenfunction is well- 
approximated by 



$o(</>) « V^cos^), 



(99) 



s=0 



i.e., the order parameter has p-wave symmetry. Con- 
tributions from higher partial waves are important only 
in the immediate vicinity of the critical angle and may 
safely be neglected otherwise (see Fig. [7]). Within the ap- 
proximation (|99p . the eigenvalue Aq equals the diagonal 
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tt/4 
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FIG. 7. (Color online) In order to solve ()97p numerically 
we expand <&o((j>) in a Fourier cosine series [the integral op- 
erator in (|97p does not couple cosine and sine series and 
it turns out that $o(<^) can be represented as a cosine se- 
ries; note that in this series expansions we retain only those 
terms that ensure the correct parity of the order parameter, 
A(£, + it) = -A(£, </>)], $oW = ^2 Eli c„ cos[(2n - 1)4>] 
with sufficiently large JV, and solve the resulting linear system 
for the coefficients c„. The above plot shows the coefficients 
ci (solid blue line), C2 (dashed green line), and C3 (dotted red 
line) as functions of 80 for do > 9 C . 



matrix element 

A «2 j gcos(0) / ^vV^f-vf')^') 
= (. 9 /7r)(4/3-3sin 2 o ). 

(100) 

In the BCS approach, therefore, for 6q > 8 C the critical 
temperature is given by 

T (BCS) _ £ (0) e -ir/g\4/3-3sin 2 e \^ 

and Eq. (I94j) determines the order parameter as 

A(e,^) = f / ^^o( P -pf)$o(A (102) 

A J 27T 

F. Critical temperature in the GM approach 

As we have pointed out earlier, within the BCS ap- 
proach it is possible to obtain only the controlling factor 
[the exponential in Eq. (|101D ] of the dependence of T c on 
g. In order to find the correct prefactor we proceed by 
substituting the ansatz 



t c = (w4 0) )t c (BCS) 



(103) 



in Eq. (j93|) . We set £ = in the resulting equation, 
multiply it by $0(0) an d take the integral over the angle 



(j). Using the normalization condition of the eigenfunction 
$0 we obtain to leading order in g 



In 



-(0) 



A o 



In 



2eT 



F 

1 ±1 



(104) 



where the correction <5Ao to the eigenvalue Ao is composed 
of four contributions, 

SX = SX { Spf} + 5\ { Sm) + £A mb) + 5\ { Q B} . (105) 

Physically, SX^ PF ^ and SX 1 ^" 1 " 1 are due to the fact that 
the pairing occurs between quasi-particles and not bare 
particles [these corrections contain the Fermi surface de- 
formation (|2"9"j) and the effective mass (J5TT) respectively, 

see discussion below], <5A mb '' has its origin in the many- 

body corrections SVd to the bare interaction, and SX 
incorporates the second order Born correction to the two- 
body vertex function. 

a. Self-energy corrections. The terms in (| 1 05|) that 
involve the self-energy T,^ and the correction to the 
chemical potential 5[i are 



SX 



(S P f) 



SX, 



(6m) 



/ 2^/ 2^ 



(0) 
P F 



dV ( P f- P f)dV (p^~pf) 



dp 



dp' 



0f 



Vo(p^-p')sW(p')A(C',0') 



f'=o 



(106) 



With the aid of Eqs. d9TJ) and (|102l) . the expression for 



SX, 



(<5?f) 



can be put in the form 



SX 



(S PF ) 



2/g*o(«/f ,*,(« 

xp-VVMpg'-pg'W*'). (107) 

where the Fermi surface deformation Spf is given by 
Eq. [the exponential smallness of T c allows us to use 
the zero temperature expression (|2"T)) for the self-energy] . 



For <5A 5m ' we find 



(6m) 



d<p' Sm((j>') 



m 

(0) (0)' 



x Vb(p^-p^)*o(A (108) 

with 5m(4>) given by Eq. pip . Performing the angular 
integrals in Eqs. (jl06[) and (|108p . we find 



SX 



(Spf) 



^ 5 2 sin 2 0o(|-isin 2 o ), 



SX 



(6m) 



f^3 2 (200 - 820sin 2 6»o + 829sin 4 (9 ). 



(109) 
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b. Many-body corrections. The leading many-body 
corrections are given by 



s \ (rob) 
oX = v 



(110) 

where the analytical expressions corresponding to the di- 
agrams ([FT) that make up SVo read 



SV (a \p,p') = [Vo(p-)} 2 

f n(q + p_/2)-n(q-p_/2) dq 
J £(q + p_/2)-£(q-p_/2)-i0(27r) s 

6V^(p,p') = -V (p-) J Vb(q- P+ /2) 

ra(q + p_/2) - n(q - p_/2) dq 
X £(q + P-/2) - £(q - P-/2) - iQ {2nf 

8Vk c \p,p') = -V (p-) J V (q + p+/2) 

ra(q + p_/2) - ra(q - p_/2) dq 
X £(q + P-/2) - C(q - P-/2) - i0 (27t)= 



, (HI) 



. (H2) 



. (H3) 



SV id) (p, p') = - J V (q - p_/2) V (q + P-/2) 
n(q + p+/2)-n(q-p+/2) dq 



£(q + p+/2) - £(q - P+/2) - *0 (2tt) 2 



. (H4) 



(we are neglecting contributions that involve excited 
states of the transverse trapping potential as they contain 
a additional factor of Ep /uio « 1). Here p± = p ± p', 

and n(p) = 6(pp —p) is the Fermi-Dirac distribution at 
zero temperature [the usage of n(p) at zero temperature 
is justified by the exponential smallness of T c ]. Perform- 
ing a numerical integration we find 



SX 



(mb) _ 2 



3^(0.37- 1.67 sin^o + 1.82 sin 4 6» ). (115) 



c. Second order Born correction. The second order 
Born correction (|85[) to the vertex function results in the 
contribution 

=/g<W / % ^H2eP,P$\pf)W) 

(116) 

to SXo. We decompose as a sum of contribu- 

tions with fixed HO quantum numbers at the inter- 
mediate state of the second order scattering process, 

^0 S) = E«,n' 5A S> where 

K : «'(Pf 0) - q) Vn,n'(q - Pf } ) rfq (n7] 
2£(q) + u; (n + n') - t0 (2tt) 2 ' 1 j 



For n = n' = we find the asymptotic expression 



SX, 



(B) ti-X) 2 



,g 2 (0.35 - 1.56sin 2 o + 1.79sin 4 6» ) 
+ g 2 01(77) (i - § sin 2 o + § sin%) - if A . (118) 



Note that despite the smallness of 77 = pf^Oj terms which 
are proportional to g 2 ln(?7) represent a small correction 
to Ao = O(g) in the limit where g/rj = Td/lo <C 1. 

We are left with the calculation of the contribution 
to the second order Born correction that involves ex- 

( B) 

cited states of the transverse trapping potential, <5A* = 

( B) ( B) 

SX — SXq . In terms of new summation indexes that 
may be interpreted as "relative" and "center of mass" 
HO quantum numbers we put it in the form 



00 N 



OA* - 2-^ 0A N+n,N-n 

N=l n=-N 

00 N 

+ E E <ttffin-l,*-». (H9) 
JV=1 n=-N+l 

i. e., we separate parts in which the sum of HO quantum 
numbers is even and odd, respectively. 

For the present purpose we rewrite the WKB matrix 
elements EU and (Ell) as 



t/ / \ / I 9 „-n 2 /2N 
V N +n,N-n{P) = \ ^T7 e ' 

V 2 TV irnv 
x [u{p) sin 2 6» - 2 P 2 (cos 6 )] 



(plo) 2 + 2N 



, (120) 



and 



» / JV+n-l,iV-n(.Pj — ( — J-J e 



x sin(2#o) «(p) 



(p(o) 2 + 2Af' 



(121) 



Inserting these expressions in Eq. (| 1 1 7[) and performing 
the integrals we find, in the limit r] — > 0, 

= -(^) 2 ^" 2/JV (1 - 4sinX + § sin%) , 

(122) 

0A N+n-l,N-n — 6 Utt JV ' e Slll^ZtfoJ . 

(123) 



In the second line we are actually restricting ourselves to 
the asymptotic behavior of S\ K N ! n _i N _ n for TV — > oo, 
which, however, gives a sufficiently accurate approxima- 
tion even for N = 1. We insert Eqs. (fl"22)) and (TT23)) in 
(fTO?)! to obtain 

* -"16x^9 

x [S t (| - 4sin 2 o + § sin 4 o ) + |5 2 sin(2# ) 2 ] , 

(124) 
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where the sums S a for a — 1 , 2 are given by 

oo 

S a = S a (N)/N 2 , (125) 

with 



N=l 



N 



Sx(N)= £ e" 



2 /N 



S 2 (N) 



n=-N 



N 

E « 

n=— W+l 



-(n-l/2) 2 /Af 



(126) 



These sums can be calculated semi-analytically with the 
result 

(5Al B) = g 2 (-0.02 - 0.01 sin 2 o + 0.03 sin 4 <9 ) . (127) 



Adding the contributions of ([TT5]) . ([TT5]l . and 

(|127l) we obtain the correction to the eigenvalue Ao, 

<5A = 9 2 (0.52 - 2.47sin 2 6» + 2.83sin 4 6» ) 

+ g 2 \n(rj) (0.25 - 1.13 sin 2 6> + 1.28 sin 4 (9 ) . (128) 



In the third term on the RHS of (|104[) we express the 
gap A(£', (f>') via (|102j) and perform a numerical integra- 
tion. We find that this term gives a contribution 



(g/X ) 2 (-0.14 + 0.63 sin 2 (9 - 0.71 sm 4 6» ) 
to ln(n7/4 0) )- 



(129) 



The combination of (fTM)) . (fT2"8"j) and (fT2"9"|) yields the 
final expression for the critical temperature in the GM 



approach 



> 



-J, 



2eTe 



(0) 



T c 
where 

f(6o) = cxp 
and 

ff(*o) 



/(0 O ) ^ o) exp(-^/ 5 |4/3 - 3 sin 2 o |) , 



0.52 - 2.47sin 2 6» + 2.83sin 4 o 



0.18- 0.81 sin 2 6» 



0.91 sin 4 6» 



0.25 - 1.13sin 2 6>o + 1.28sin 4 6> 
0.18 - 0.81 sin 2 6> + 0.91 sin 4 (9 ' 



(130) 
(131) 

(132) 



Fig. [5] shows the critical temperature as a function of 
the tilting angle for values of g and r\ that correspond 
to a gas of polar KRb molecules. The corresponding 
eigenfunction $o(0) for the order parameter is given by 

HMD- 



IX. SUMMARY AND CONCLUSION 

In this paper, we studied various properties of a quasi- 
2D dipolar Fermi gas. We found that the normal phase 
is characterized by an anisotropic Fermi surface, which - 
in order to minimize the interaction energy - is elongated 



2 - 



1 - 



tt/4 



tt/2 



FIG. 8. The critical temperature as a function of 0q for g ~ 
1.2 and r\ w 0.2. These values correspond to a gas of KRb 



molecules, with U2d = 4- 
discussion in Sec. IIXI 

tt/2 



3tt/ 



10 8 cm"- 1 and lj = 2tt ■ 100 kHz, see 
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FIG. 9. (Color online) Phase diagram of the quasi-2D dipo- 
lar Fermi gas at T = 0. For < g < 1.45 and small tilting 
angles the system is a normal Fermi liquid (NFL). The tran- 
sition to the superfluid phase (SF) occurs at the critical angle 
8 C sa 0.7217. In our perturbative approach, which is - strictly 
speaking - valid only for g <C 1, this value is independent of 
g. At moderately strong interactions there appear a density 
wave phase (DW) in the lower part of the phase diagram and 
a long wavelength instability (LWI) in the upper part. 



along the projection of the dipoles on the plane of con- 
finement (here: the x-axis). Consequently, the dispersion 
relation of single-particle excitations is anisotropic too, 
with an increased effective mass in the x-direction and a 
decreased effective mass in the y-direction (see Fig. [5]). 

Long wavelength collective excitations (zero sound), 
corresponding to deformations of the equilibrium Fermi 
surface, can propagate through the medium only when 
the tilting angle exceeds a value of 9 ~ 0.5, which de- 
pends only weakly on g for (at least) 0.1 < g < 2. In 
this regime we found two distinct modes, symmetric and 
antisymmetric about the propagation direction, respec- 
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tively. 

The experimental observation of zero sound requires T 
so low that the attenuation is negligible. Note, however, 
that for 6*o > 6* c ~ 0.7 (the critical angle for the superfluid 
transition) we also require T > T c for the observation 
of zero sound, although under certain conditions, zero 
sound can still persist below T c (see discussion in [33[ 
and references given therein). 

At higher values of g and for large tilting angles, we 
found complex eigenfrequencies of the long wavelength 
collective modes indicating the instability of the system 
towards collapse in this part of the phase diagram, see 
Fig. O Instabilities with finite momentum that drive the 
system towards a state with periodic density modulation, 
which breaks translational symmetry, appear at small 
tilting angles. In the isotropic case with 9o = (dipolcs 
are aligned perpendicular to the plane of confinement) 
also rotational symmetry is broken. This phase can be 
detected experimentally using Bragg spectroscopy. 

Finally, let us discuss prospects for observing the su- 
perfluid transitions in the quasi-2D dipolar Fermi gas. In 
current experiments with KRb at JILA, electric dipole 
moments of d ~ 0.3 D are available. Then, assuming 
an area density of n-io ~ 4 • 10 8 cm 2 corresponding to 
Ep ~ 100 nK, the dimensionless coupling constant takes 
the value g « 1.2. Note that strictly speaking this 
value of g lies beyond the weak coupling regime in which 
Eq. (I130[) is valid. However, this formula still provides an 
estimate for the onset of superfluidity in the intermediate 
coupling regime. 

The strong confinement condition is realized for a 
transverse trapping frequency of luq = 2n ■ 100 kHz, so 
that r] w 0.2. For the above value of g the system is sta- 
ble up to a tilting angle of 8q 1.3 (see Fig. [!)]). Then we 
have T c ss 0.04 « 4nK, which is comparable to the 
values of T r that are to be expected in bilayer systems 
(seeRefs. [lMl EMI). 
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Appendix A: Calculation of the HF self-energy 

In order to solve Eqs. (|2"2"]l and ([2~4"]l for pp it is conve- 
nient introduce another unknown function /(</>) accord- 
ing to 



The expansion of / in a Fourier series takes the form 



f{<t>)=<p F {<t>f/pf 2 



f(&) = }] c„cos(2n0), 



n£ 



(A2) 



i.e., the constant term is absent, as is immediately ap- 
parent by inserting the definition (|Aip in Eq. (1231) . More- 
over, due to the symmetry of the problem, terms which 
are proportional to sin(n</>) or cos(n0) with odd n do not 
appear. 

We solve for the Fourier coefficients c„ by iteration: As 
initial values we take the results from first order pertur- 
bation theory (see Sec. W} . 



„(D = ) T^9 sin O for n = 1, 

|0 for n = 2,3, . 



(A3) 



The iteration scheme consists in repeatedly carrying 
out the integrals [which result from the combination of 
Eqs. (HU), (El, HID), and (TA2l) ] 



1 



dx I d(j)d<f>' cos(2n0)a; l + / (l) (0') 



x v V (pf ^1 + /«(«£) p - pfxyjl + 



IP' 



(A4) 

for given values of g and 9q, where we keep coefficients up 
to n = 4 [i. e., terms in the series (|A2[) up to cos(80)] and 

/^(</>) is obtained by replacing c„ in Eq. (|A2[) by c£ ■ 
A measure for the convergence of this procedure is the 
relative change in / in one step of the iteration process, 



[ /(*+!) - /(*)| 



\ 



E, 



(i+i) _ (») 

Cn 



J2 n Cn 



(A5) 



and we terminate the iteration when this quantity drops 
below 10~ 3 . 

From the final result for / we immediately obtain pp by 
inverting Eq. (IA1|) . Then, S and e follow from Eq. (Til)!) , 
and by taking the gradient in (1191) we obtain the Fermi 
velocity (l2l)f as 



1 

V F = — 

m 



r dS' f p ' F 
VF J H Q *VfW (PF-p') 



(A6) 



Appendix B: Numerical solution to the 
Bethe-Salpeter equation 



To solve Eq. (|32|) numerically it is advantageous to 
rewrite it in a form that is symmetrized with respect to 
the transferred momentum q. To this end, we introduce 
the shifted function Xq(p) = x(P — <l/2;q). Then we 
have 



f dp' ~ . 

Xq(P) = J 7^)2 r P h (P - P >q) 



1. 



(Al) 



{2-Kf 

n(p' - q/2) - n(p> + q/2) 
w + e(p'-q/2)-e(p' + q/2) 



Xq(p')- (Bl) 
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The factor n(p' — q/2) — n(p' + q/2) in the numerator 
on the RHS restricts the area of integration as depicted 
in Fig. [TDJ Rewriting the integral in polar coordinates 
p' = p' (coscf)', sin^'), and in terms of a new function 

" q(P) S c + s(p-q/2) P - £ (p + q/2)' (B2) 



FIG. 11. (Color online) Top: Solid and dashed lines corre- 
spond to Vn%(p) with n — 0,2, ... ,10, calculated exactly and 
in WKB approximation, respectively. (We do not use the 
WKB approximation in the case n = 0, hence there is no 
dashed line corresponding to uppermost solid line.) Bottom: 

Same as above, for V^ I } l , 1 (p) with n = 0, 2, . . . , 10. In both 
figures higher values of n correspond to darker colors. 



Appendix C: The matrix elements V n ^ n > 



Eq. (IBip becomes (for the definition of the angles (j>± that 
limit the (//-integration see Fig. [TU] p± (</>') are the angle- 
dependent upper and lower boundaries for the integration 
over p') 

[u + e(p - q/2) - e(p + q/2)] i/ q (p) 

i d6' r p +^"> 

= -k t- d P 'p'r ph (p-p', q )^(p') 

2ir J (j> _ Zir J p _(p) 
i r^+<l>+ aaj rP+(4>') 

-7T it dp'p'T ph (p-p',<l)^(p')- 

(B3) 

We introduce a new variable x € [0, 1] that parametrizes 
the momentum as p(x, 4>) — P-{4>) + \p+(4>) ~ P-(4>)] x , 
and discretize the resulting integrals over x' and <j>' ac- 
cording to the trapezoidal quadrature rule in two dimen- 
sions. Here we choose a number of 20 grid points in the 
variable x 1 and 140 grid points in <j)' . Then, for the val- 
ues of the parameters of Fig. H2 an addition of 10 grid 
points in either variable leads to an absolute change in 
the result for w/ep^ that is less than 10~ 2 . 



In this appendix we present explicit expressions for the 
matrix elements ([5]) with — — 0, which appear in 
the second order term (|85|) in the Born series for the 
two-body vertex function. Omitting the momentum in- 
dependent contribution (see footnote [IJ , with the aid of 
the integral table [5l| we find, for n + n' e 2No, 

V n y{p) = (2"+"W!)-V2 1X1+2^) 

x (d 2 /l ) [u{p)sin 2 8 -2P 2 (cos9 )] 

x ( p l ) 1 + n + n ' e P 2l o /2 T(^^,p 2 l 2 /2), (CI) 

where u(p) = (p 2 —p 2 )/p 2 , and T(z) and T{y,z) are the 
complete and incomplete Gamma functions respectively. 
For n + n' € 2No + 1, with v(p) = pi/p, we have 

V n ,„> (p) = [(-i)»-«'/2 n + n ' n \ n'!] 1 / 2 r(l + 2±ai) 
x (2d 2 //o)w(p)sin(26»o) 
x {pl ) 1+nW e p2l ° /2 T{-^,p 2 ll/2). (C2) 

The relatively complicated functional dependence of 
the matrix elements V n ^ n ' on the momentum p prevents 
us from calculating the integrals in (111 7[) analytically. We 
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avoid this problem by using WKB matrix elements [52J 
instead of the exact ones. 

The WKB eigenfunctions of the one-dimensional HO, 
which give good approximations to the exact eigenfunc- 
tions for n > 1, are Q 



4> n {z) = {-l) n ^2mu /irpn(z) cos[* n (z) - tt/4], (C3) 

where the position dependent momentum p n {z) and the 
phase 9 n (z) are given by 



Pn (z) = V2m[£^- 
#n(z) = / dz'p n (z'), 



(C4) 
(C5) 



with the potential V(z) = muS^z 1 j2 and the HO energy 
levels £„ = ^o(" + 1/2)- Above expression for 0„(z) is 
valid in the classically allowed region [i. e., for values of z 
such that E n > V(z)] and far away from the two turning 
points at ±z„ = ±loy/2n + 1, at which the kinetic energy 
is zero, p n (±z n ) = , or, equivalently, E n — V(±z n ). 

The HO ground state wave function decays exponen- 
tially on a scale that is set by lo. Therefore, the integra- 
tion in the matrix elements V n ,n' is essentially restricted 
to the interval \z\ < Iq. For these values of z we may use 
the approximations 

Pn(z) ~ p n (0) = p n , ^ n(z) ~ p„Z + (%/A){2n + 1). 

(C6) 

The WKB wave function (IC3I) then becomes 



Consequently, for the matrix elements V n>n > in WKB ap- 
proximation we find 

V n , n ,(p) « 2(n + n')- 1/2 e-("-"') 2 / 4 ("+"')(d 2 // ) 

x [ U (p)sin 2 o - 2P 2 (cos0 o )] {ph yf n + nl (C8) 

for even n + n', and 



K,n'(p) ~ 4i n_n ' e -("-"') 2 /4(i+«+«') 
/n 



X^sm(2go)^(p) 73 ^ , (C9) 



{pl ) 2 + 1 + n + n' 



for odd n + nf. 

In order to conveniently compare the exact results for 
V n , n > with those obtained in WKB approximation we in- 
troduce reduced matrix elements which are real, dimen- 
sionless and depend only on the magnitude p of the mo- 
mentum. For even n + ml the reduced matrix element 
Van' (p) i s defined via the relation 



V n ,n>(p) = £ [u(p)sm 2 9 - 2P 2 (cos6 )] V£l(p), 



and, for odd n + n' , we set 



(CIO) 



(Cll) 



4>n{z) w (-l) n x/2muj /irp n cos(p n z + nir/2). (C7) 



Figure [TT] compares the reduced matrix elements in WKB 
approximation with the exact ones. Quantitative agree- 
ment improves with increasing HO quantum numbers 
(corresponding to darker colors in Fig. 1111) and is, how- 
ever, satisfactory even for n, n' > 1. 
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